###
### Code to replicate tables and figure in the main text
###




###
### Tables 1 and 3
###
library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=3)

# multi-level negative binomial model
nb.model <- "
data{
  int N;
  int K;

  int y[N];
  
  matrix[N,K] X;

  int<lower=0> P;
  int<lower=0, upper=P> party[N];

  int<lower=0> S;
  int<lower=0, upper=S> province[N];
}

parameters{
  real alpha;
  vector[K] beta;

  vector[P] nu_p;
  vector[S] nu_s;
  
  real<lower=0> phi;
  
  real<lower=0, upper=100> sigma_p;
  real<lower=0, upper=100> sigma_s;
}

model{
  y ~ neg_binomial_2_log(alpha + X*beta + nu_p[party] + nu_s[province], phi);

  nu_p ~ normal(0, sigma_p);
  nu_s ~ normal(0, sigma_s);
}
"

# multi-level linear model
normal.model <- "
data{
  int N;
  int K;

  vector[N] y;

  matrix[N,K] X;

  int<lower=0> P;
  int<lower=0, upper=P> party[N];

  int<lower=0> S;
  int<lower=0, upper=S> province[N];
}

parameters{
  real alpha;
  vector[K] beta;

  vector[P] nu_p;
  vector[S] nu_s;

  real<lower=0, upper=100> sigma_y;

  real<lower=0, upper=100> sigma_p;
  real<lower=0, upper=100> sigma_s;
}

model{
  y ~ normal(alpha + X*beta + nu_p[party] + nu_s[province], sigma_y);

  nu_p ~ normal(0, sigma_p);
  nu_s ~ normal(0, sigma_s);
}
"

sm.nb <- stan_model(model_code=nb.model)
sm.normal <- stan_model(model_code=normal.model)

# load data
load("NationalLegislatorsData.RData")
colnames(df)

# dummy for governing party
df$gov.party <- 0
df$gov.party[df$party=="PML(N)"] <- 1
table(df$gov.party)

# party and province IDs
df$partyID <- as.numeric(as.factor(df$party))
df$provinceID <- as.numeric(as.factor(df$province))

# three types of legislators 
# type2 = SMDP female
# type3 = RS female
# type4 = RS religious minority
df$type2 <- 0
df$type2[df$type==2] <- 1
df$type3 <- 0
df$type3[df$type==3] <- 1
df$type4 <- 0
df$type4[df$type==4] <- 1

round(cor(df[,c("type2", "type3", "type4")]), 2)

prop.dat1 <- list(N=nrow(df),
                  y=df$n.proposal, 
                  X=cbind(df$reserve, df$term, df$committee_chair),
                  K=3,
                  P=length(unique(df$partyID)),
                  party=df$partyID,
                  S=length(unique(df$provinceID)),
                  province=df$provinceID)

# table 1, model 1
prop.fit1 <- sampling(sm.nb, data=prop.dat1, iter=15000, warmup=10000, chains=3,
                      seed=1234)

prop.dat2 <- list(N=nrow(df),
                  y=df$n.proposal, 
                  X=cbind(df$type4, df$type3, df$type2, 
                          df$term, df$committee_chair),
                  K=5,
                  P=length(unique(df$partyID)),
                  party=df$partyID,
                  S=length(unique(df$provinceID)),
                  province=df$provinceID)

# table 1, model 2
prop.fit2 <- sampling(sm.nb, data=prop.dat2, iter=15000, warmup=10000, chains=3,
                      seed=1234)

cos.dat1 <- list(N=nrow(df),
                 y=df$n.cosponsor, 
                 X=cbind(df$reserve, df$term, df$committee_chair),
                 K=3,
                 P=length(unique(df$partyID)),
                 party=df$partyID,
                 S=length(unique(df$provinceID)),
                 province=df$provinceID)

# table 1, model 3
cos.fit1 <- sampling(sm.nb, data=cos.dat1, iter=15000, warmup=10000, chains=3,
                     seed=1234)

cos.dat2 <- list(N=nrow(df),
                 y=df$n.cosponsor, 
                 X=cbind(df$type4, df$type3, df$type2, 
                         df$term, df$committee_chair),
                 K=5,
                 P=length(unique(df$partyID)),
                 party=df$partyID,
                 S=length(unique(df$provinceID)),
                 province=df$provinceID)

# table 1, model 4
cos.fit2 <- sampling(sm.nb, data=cos.dat2, iter=15000, warmup=10000, chains=3,
                     seed=1234)

bet.dat1 <- list(N=nrow(df),
                 y=log(df$between+1), 
                 X=cbind(df$reserve, df$term, df$committee_chair),
                 K=3,
                 P=length(unique(df$partyID)),
                 party=df$partyID,
                 S=length(unique(df$provinceID)),
                 province=df$provinceID)

# table 3, model 1
bet.fit1 <- sampling(sm.normal, data=bet.dat1, iter=15000, warmup=10000, chains=3,
                     seed=1234)

bet.dat2 <- list(N=nrow(df),
                 y=log(df$between+1), 
                 X=cbind(df$type4, df$type3, df$type2, 
                         df$term, df$committee_chair),
                 K=5,
                 P=length(unique(df$partyID)),
                 party=df$partyID,
                 S=length(unique(df$provinceID)),
                 province=df$provinceID)

# table 3, model 2
bet.fit2 <- sampling(sm.normal, data=bet.dat2, iter=15000, warmup=10000, chains=3,
                     seed=1234)

ego.dat1 <- list(N=nrow(df),
                 y=log(df$ego_frac+1), 
                 X=cbind(df$reserve, df$term, df$committee_chair),
                 K=3,
                 P=length(unique(df$partyID)),
                 party=df$partyID,
                 S=length(unique(df$provinceID)),
                 province=df$provinceID)

# table 3, model 3
ego.fit1 <- sampling(sm.normal, data=ego.dat1, iter=15000, warmup=10000, chains=3,
                     seed=1234)

ego.dat2 <- list(N=nrow(df),
                 y=log(df$ego_frac+1), 
                 X=cbind(df$type4, df$type3, df$type2, 
                         df$term, df$committee_chair),
                 K=5,
                 P=length(unique(df$partyID)),
                 party=df$partyID,
                 S=length(unique(df$provinceID)),
                 province=df$provinceID)

# table 3, model 4
ego.fit2 <- sampling(sm.normal, data=ego.dat2, iter=15000, warmup=10000, chains=3,
                     seed=1234)

#save(prop.fit1, prop.fit2, cos.fit1, cos.fit2,
#     bet.fit1, bet.fit2, ego.fit1, ego.fit2, file="PakistanStan.RData")

# summary
load("PakistanStan.RData")

# table 1, model 1
round(summary(prop.fit1)$summary[c(2:4,31,32,30),c(1,3)], 3)

# table 1, model 2
round(summary(prop.fit2)$summary[c(2:6,33,34,32),c(1,3)], 3)

# table 1, model 3
round(summary(cos.fit1)$summary[c(2:4,31,32,30),c(1,3)], 3)

# table 1, model 4
round(summary(cos.fit2)$summary[c(2:6,33,34,32),c(1,3)], 3)

# table 3, model 1
round(summary(bet.fit1)$summary[c(2:4,31,32),c(1,3)], 3)

# table 3, model 2
round(summary(bet.fit2)$summary[c(2:6,33,34),c(1,3)], 3)

# table 3, model 3
round(summary(ego.fit1)$summary[c(2:4,31,32),c(1,3)], 3)

# table 3, model 4
round(summary(ego.fit2)$summary[c(2:6,33,34),c(1,3)], 3)




###
### Table 2
###
library(statnet);library(latentnet)

load("NationalLegislatorsData.RData")
load("CosponsorshipNetwork.RData")

# dummy for governing party
df$gov.party <- 0
df$gov.party[df$party=="PML(N)"] <- 1
table(df$gov.party)

cos.c <- as.network(cosponsor, loops=FALSE, directed=FALSE,
                    ignore.eval=FALSE, names.eval="count")
set.vertex.attribute(cos.c, "reserved", df$reserve+1)
set.vertex.attribute(cos.c, "type", df$type)
set.vertex.attribute(cos.c, "term", df$term)
set.vertex.attribute(cos.c, "committee_chair", df$committee_chair)
set.vertex.attribute(cos.c, "gov.party", df$gov.party)
set.vertex.attribute(cos.c, "party", df$party)
set.vertex.attribute(cos.c, "province", df$province)

lsm.f3 <- ergmm(cos.c ~ nodefactor("reserved", base=1) + nodefactor("gov.party")
                + nodecov("term") + nodefactor("committee_chair") 
                + nodematch("party") + nodematch("province")
                + nodemix("reserved", base=c(1,1))
                + euclidean(d=3),
                family="Bernoulli", response="count", tofit="mle")

load("LatentPredictions.RData")

# table 2
summary(lsm.f3)




###
### Figure 1
###
library(igraph)

load("NationalLegislatorsData.RData")
load("CosponsorshipNetwork.RData")

# remove isolates
num <- which(rowSums(cosponsor) != 0)
cosponsor <- cosponsor[num, num]

net <- graph.adjacency(cosponsor, mode="undirected", diag=FALSE)
V(net)$color <- ifelse(df$reserve[num]==1, "black", "white")
V(net)$label <- NA
V(net)$size <- ifelse(df$between[num] < 1000, 4, 8)
E(net)$width <- 0.4
E(net)$color <- "gray70"

# figure 1
#pdf("pak_net.pdf")
#tiff("pak_net.tif", width=2000, height=2000, res=400)

set.seed(1234)
par(mar=c(1,1,1,1))
plot(net)
legend("bottomright", c("RS Legislator", "Non-RS Legislator"),
       col=c("black", "black"), pch=c(19,21), cex=1.1)

#dev.off()